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ABSTRACT 

A  detailed  parametric  study  of  the  variables  critical  to  a  numerical 
model  of  the  interaction  of  a  free  surface  with  an  ascending  vortex 
pair  has  been  carried  out.  The  model  is  based  on  a  discrete  vortex 
representation  of  the  free  surface  and  the  use  of  the  appropriate 
boundary  conditions.  A  study  was  performed  to  minimize  the  numeri- 
cal Helmholtz  instabilities  inherent  in  the  numerical  model  for  lower 
Froude  numbers.  A  method  of  free-surface  re-discretization  was 
developed  and  implemented.  Additionally,  an  investigation  of  the 
motion  of  a  vortex  pair  approaching  and  passing  through  a  sharp 
salt/freshwater  density  interface  was  conducted.  Experiments  were 
performed  in  a  large  basin  and  the  vortices  were  generated  with  two 
counter-rotating  plates.  The  rise  and  intrusion  of  the  Kelvin  oval  into 
the  upper  layer  of  lower  density  fluid  was  recorded  on  a  video  tape  and 
then  analyzed  through  the  use  of  a  Motion  Analysis  system. 
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I.    INTRODUCTION 

The  understanding  of  the  temporal  development  of  the  surface 
signatures  resulting  from  the  interaction  of  an  ascending  inclined  vor- 
tex pair  with  the  free  surface  is  of  special  importance.  When  a  trailing 
vortex  pair,  generated  by  a  lifting  surface,  approaches  the  free  surface, 
with  or  without  mutual  induction  instability  and /or  vortex  breakdown, 
the  vortices  and /or  the  crude  vortex  rings  give  rise  to  surface  distur- 
bances (scars  and  striations).  The  striations  are  essentially  three- 
dimensional  free-surface  disturbances,  normal  to  the  direction  of  the 
motion  of  the  lifting  surface,  and  come  into  existence  when  the  vortex 
pair  is  at  a  distance  equal  to  about  one  initial  vortex  separation  from 
the  free  surface.  The  scars  are  small  free-surface  depressions  on 
either  side  of  a  nearly  symmetrical  hump,  directly  above  the  plane  of 
symmetry  of  the  vortex  pair.  They  come  into  existence  towards  the 
end  of  the  pure  striation  phase  and  when  the  vortices  are  at  a  distance 
equal  to  about  60  percent  of  the  initial  vortex  separation  from  the  free 
surface.  The  scars  are  created  and  driven  by  the  vortices  and  both 
move  outward.  Subsequently,  the  vortices  diffuse  rapidly  and  their 
circulation  decreases  slowly  at  first  (prior  to  the  occurrence  of  scars) 
and  rapidly  thereafter  due  to  the  counter-sign  vorticity  generated  in 
the  overlapping  regions  of  the  vortex  pair  in  the  recirculation  cell  and 
at  the  free  surface. 


The  surface  disturbances  or  the  local  roughening  of  the  free  sur- 
face, resulting  from  the  vortex  motion  produced  by  submerged  foils  or 
counter-rotating  propellers,  are  visible  because  the  particle  motions 
they  create  as  they  propagate  produce  a  straining  field  at  the  surface 
which  modifies  the  physical  characteristics  of  the  sea  surface 
(reflection,  refraction,  and  the  scattering  of  sound,  heat,  and  light).  In 
fact,  in  many  photographs  of  the  sea  surface  taken  from  aircraft  or 
earth  resources  technology  satellites,  using  passive  or  active  sensors, 
the  differential  roughening  of  the  free  surface  can  be  seen  clearly. 

Previous  studies  of  a  vortex  (or  vortex  pair)  in  ground-  or  free- 
surface  effect  have  been  generally  concerned  with  the  boundary  layer 
development,  generation  of  oppositely-signed  vorticity,  the  flow 
resulting  from  a  single  vortex  held  fixed  in  a  uniform  flow  and  the  so- 
called  "rebounding"  phenomenon,  according  to  which  the  vortices 
move  at  first  toward  the  boundary,  as  on  an  inviscid  trajectory,  and 
then  away  from  it.  The  boundary  surface  may  be  either  a  rigid  bound- 
ary at  which  the  no-slip  condition  is  satisfied  or  a  free  surface  corre- 
sponding to  zero  shear  stress.  In  either  case,  the  boundary  surface 
(sometimes  erroneously  referred  to  as  a  vorticity  sink)  gives  rise  to 
oppositely  signed  vorticity  and  the  net  circulation  of  the  vortices 
decreases  with  time. 

The  interaction  of  vortices  with  a  fluid  interface  has  been  under 
investigation  at  the  Naval  Postgraduate  School  by  Sarpkaya  and  his  stu- 
dents for  the  past  several  years  [Refs.   1-11].    The  emphasis  of  these 


investigations  has  been  the  identification  of  the  specific  mechanisms 
which  lead  to  the  demise  of  the  vortices  and  the  prediction,  through 
experiments  and  analysis,  of  the  motion  of  vortices  and  the  form  of  the 
free  surface. 

Sarpkaya  and  Henderson's  [Ref.  9]  theoretical  analysis  of  the  scars 
created  by  the  trailing  vortices  was  carried  out  assuming  the  vortices 
to  be  two-dimensional  and  the  free  surface  to  be  a  rigid  plane.  For 
small  Froude  numbers  (Fr  =  V0/ygb0) ,  (where  V0  is  the  initial 
mutual  induction  velocity  of  the  vortices  and  b0  is  the  initial  vortex 
separation),  each  scar  front  was  shown  to  coincide  with  the  stagnation 
point  on  the  Kelvin  oval,  formed  by  one  of  the  pair  of  the  trailing 
vortices  and  its  image.  For  Froude  numbers  larger  than  about  0.15, 
not  only  the  deformation  of  the  free  surface  but  also  the  nonlinear 
interaction  between  the  said  deformation  and  the  motion  of  the 
vortices  is  significant.  Thus,  the  free  surface  may  no  longer  be 
assumed  rigid.  However,  the  analysis  of  a  three-dimensional 
interaction  between  the  surface  and  the  inclined  vortex  pair  is 
extremely  difficult.  It  is  in  view  of  this  fact  that  it  was  decided  to 
generate  nearly  two-dimensional  vortex  pairs  at  relatively  high  Froude 
numbers  and  to  carry  out  a  numerical  analysis  assuming  a  deformable 
free  surface.  Most  recently,  Leeker  [Ref.  Ill  continued  the  work  on 
the  interaction  of  a  vortex  pair  with  the  free  surface.  His  task  was  to 
improve  upon  the  discrete  vortex  model,  originally  constructed  by 
Elnitsky  [Ref  10],  by  reducing  or  eliminating  the  numerical 
instabilities  present  at  lower  Froude  numbers.     This  improvement 


would  enable  the  model  to  more  closely  represent  the  results 
observed  in  the  laboratory  by  allowing  the  calculations  to  progress 
further  in  time. 


The  present  investigation  has  a  two-fold  purpose: 

1.  To  conduct  a  detailed  parametric  study  of  the  variables  affecting 
the  performance  of  the  discrete  vortex  model  and  identify  the 
variables  (if  any)  that  adversely  or  positively  affect  the  perfor- 
mance of  the  model  with  particular  emphasis  on  their  effects  on 
observed  numerical  instabilities.  Once  it  has  been  determined 
how  sensitive  the  numerical  model  is  to  a  variety  of  variables, 
efforts  are  made  to  attempt  a  qualitative  optimization  of  these 
variables  under  specific  conditions.  Armed  with  the  information 
gained  in  the  parametric  study,  additional  efforts  are  made  to 
obtain  methods  that  remove  or  neutralize  numerical  instabilities 
present  at  lower  Froude  numbers.  Several  smoothing  techniques 
and  re-discretization  schemes  are  attempted. 

2.  To  investigate  the  motion  of  a  vortex  pair  approaching  and 
intruding  through  a  sharp  density  interface.  Utilizing  an 
improved  test  apparatus  and  a  sharp  salt/fresh  water  interface,  a 
complete  investigation  was  conducted  to  analyze  the  interaction 
of  the  vortices  and  the  Kelvin  oval  with  the  interface.  The  motion 
of  the  vortex  centers  was  recorded  on  video  tape  and  analyzed 
through  the  use  of  a  motion  analysis  system. 


H.   PARAMETRIC  STUDY  OF  THE  NUMERICAL  MODEL 

A      INTRODUCTION 

The  discrete  vortex  model  used  in  this  investigation  and  in  the 
previous  work  is  discussed  by  Elnitsky  and  Leeker  in  detail  [Refs.  10, 
11].  In  summary,  discrete  vortices  are  placed  along  the  free  surface  to 
form  a  vortex  sheet.  The  inter-vortex  spacing  may  vary  as  desired, 
emanating  from  a  pre-selected  center  of  vortex  distribution  along  the 
free  surface.  After  the  program  is  executed,  each  individual  surface 
vortex  is  convected  utilizing  the  velocities  imparted  to  it  by  all  other 
vortices  and  the  main  rising  vortex  pair.  Execution  continues  until 
numerical  Helmholtz  instabilities  develop  and  chaos  ensues. 

The  complex  function  representing  all  of  the  vortices  in  the  sys- 
tem is  given  by: 

ln(z  -zo)      ln(z  +  z0)      ^rkln(z-zk)      ^iy^z+zj 
w  =  ~ — oT^ — " oT^ ^ oT^ +  ^ oT^ (1) 

ZlK  ZlK  ZlTZ  ,  Z17U 

k=l  k=l 

The  strengths  of  the  free-surface  vortices  are  normalized  by  the 
strength  of  the  main  vortex.  Spatial  units  are  normalized  by  the  dis- 
tance b0,  the  initial  spacing  between  the  two  main  vortex  centers. 
Velocities  are  normalized  by  V0  =T0/2b0  and  time  is  normalized  by 
T*=V0t/b0. 

The  normalized  boundary  condition  at  the  free  surface  is  then 
given  by: 


DT*      2        Fr2  l    J 

and  the  induced  velocity  becomes: 

dw  _     _. i i y      1X  m         y        1X  m  ^) 

dt~  ~2n(z-z0)      2k(z  +  z0)      ^2tc(z-Zi)     ^2k{z  +  z) 

There  are  other  numerical  methods  available  to  solve  the  problem 
of  the  interaction  of  a  free  surface  with  a  vortex  pair,  and  Elnitsky  dis- 
cusses them  in  great  detail  [Ref  10].  The  principal  advantages  of  the 
Boundary  Integral  Equation  method  utilized  in  this  and  preceding 
investigations  are  the  ability  to  discretize  the  free  surface  without  a 
loss  of  accuracy  and  the  ability  to  solve  the  boundary  conditions 
exactly.  Another  advantage  that  this  method  enjoys  is  the  relative 
computer  efficiency  with  which  it  performs  with  respect  to  total  com- 
putational time.  However,  in  order  to  utilize  this  model  to  accurately 
predict  the  behavior  of  the  motion  of  a  vortex  pair  in  nature,  it  is  nec- 
essary that  a  means  be  discovered  to  minimize  the  impact  of  the 
numerical  instabilities.  Finite  Element  and  Finite  Difference  tech- 
niques were  also  considered  as  a  means  of  solving  this  problem  but 
were  determined  to  be  too  cumbersome  and  computationally 
intensive. 

The  initial  step  in  the  task  of  minimizing  numerical  instabilities 
within  this  distributed  vortex  model  is  to  identify  which  variables  most 
affect  the  performance  of  the  numerical  code.  In  concert  with  this 
task  is  the  requirement  to  qualitatively  assess  the  relative  sensitivity  of 


the  code  to  each  changing  variable.  Finally,  after  each  variable  has 
been  qualitatively  optimized,  additional  methods  must  be  imple- 
mented to  further  improve  the  numerical  stability  of  the  model  for 
lower  Froude  numbers.  The  effect  of  the  following  parameters  were 
investigated: 

5:  desingularization  parameter 

N:  total  number  of  free  surface  discrete  vortices 

xmax:      maximum  x/b0  position 

xv:  center  of  the  surface  vortex  distribution  (x/b0) 

Upon  the  completion  of  the  parametric  study,  a  re-discretization 
method  was    devised  and  implemented. 

R      DESINGULARIZATION  PARAMETER  (5) 
1.      Introduction 

Leeker  [Ref  11]  discussed  a  desingularization  method  intro- 
duced by  Rosenhead  [Ref  12]  which  uses  a  smoothing  parameter  to 
remove  the  singularity  at  the  center  of  a  given  vortex  and  to  minimize 
the  induced  velocity  of  a  vortex  on  its  nearest  neighbors.  Utilizing  this 
method  of  desingularization,  equation  (3)  becomes: 

dw  m        =  Kz-zo)2  Kz  +  2o)2 

dt      U   1V        2  7t(z-zo)[(z-zo)2  +  52]'f  2K(z  +  z0i(z+z0)2+52] 

+  £  ^(z-z,)2  _£  iiyz+z/  (4) 

k=i    2tc(z-z,)[(z-z1)2  +  5  ]    k=i    2tc(z  +z,)[(z  +  z1)    +o\ 


Clearly,  this  method  of  desingularization  eliminates  the 
unnatural  infinite  velocity  at  the  vortex  center  and  allows  the  operator 
to  more  closely  space  individual  vortices  at  the  free  surface.  Reduced 
inter-vortex  spacing  minimizes  the  leakage  problem  arising  from 
widely  spaced  vortices,  yet  without  the  desingularization  parameter, 
this  close  spacing  would  immediately  lead  to  individual  surface  vor- 
tices rotating  about  one  another,  creating  irreparable  numerical  insta- 
bilities. It  is  clear  that  as  the  spacing  between  vortex  elements  is 
increased,  the  effect  of  8  on  the  performance  of  the  code  would 
diminish  accordingly.  In  order  to  more  fully  understand  the  impact  of 
this  de-singularization,  a  systematic  parametric  study  was  undertaken. 
From  earlier  investigations  by  Elnitsky  and  Leeker  [Refs.  10, 
11]  the  optimum  parameters  for  executing  the  numerical  model  and 
ensuring  that  it  would  progress  an  adequate  amount  of  time  were 
qualitatively  determined  to  be: 

Number  of  surface  vortices  (N)  76 

Froade  number  (Fr)  .35 

Starting  depth  (d0)  -3.0 

Starting  time  (T*)  -3.00 

2.     Application 

The  smoothing  parameter  8  was  varied  from  0.05  to  0.15  and 
the  program  was  executed  from  a  non-dimensional  starting  time  of 
-3.00  The  choice  of  8  was  governed  by  the  number  of  surface  ele- 
ments that  would  be  affected  at  one  time  by  a  given  vortex.  If  8  is  cho- 
sen too  large  with  respect  to  the  quantity  (z-Zi),   the  program  will 
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essentially  ignore  several  consecutive  vortex  contributions  in  its  cal- 
culation of  the  velocity  field;  this  is  particularly  prevalent  at  the  early 
stages  of  the  motion  where  the  vortex  spacing  is  comparatively  small 
near  the  axis  of  symmetry.  Ideally,  the  goal  was  to  avoid  having  5  mask 
more  than  two  consecutive  surface  vortices  at  the  beginning  of  the 
distribution  and  yet  purposely  mask  vortices  where  the  numerical 
instabilities  form.  This  allows  the  calculation  to  proceed  with  each 
vortex  having  only  a  minor  effect  on  its  nearest  neighbors  but  still 
having  the  full  effect  on  more  distant  vortices.  As  a  result,  the  opera- 
tor may  space  vortices  very  close  together  and  yet  still  minimize  the 
instabilities  which  would  otherwise  occur. 

Table  1  summarizes  a  representative  number  of  executions  of 
the  numerical  model  while  holding  all  other  variables  constant  and 
varying  6  from  .05  to  .15.  The  breakdown  time  represents  the  last 
reliable  time  (±  five  time  steps)  where  the  free  surface  is  not  unnatu- 
rally formed.  As  can  be  seen  clearly  from  Figure  1,  the  lowest  value  of 
6  (.05)  is  unacceptable.  The  free  surface  deforms  in  an  entirely  unnat- 
ural fashion  at  all  non-dimensional  times  later  than  -1.4.  Run  2 
demonstrates  significant  improvement  with  only  a  very  minor  increase 
in  6  (Figure  2).  For  subsequent  runs  with  a  slightly  increasing  6  value, 
the  gain  in  extended  execution  time  is  minimized  and  appears  to  con- 
verge somewhat  to  a  non-dimensional  time  of  about  -1.0.  Figures  3 
through  7  illustrate  the  form  of  the  free  surface  for  each  run  at  the  last 
reliable  time  just  prior  to  the  halt  of  the  execution  due  to  numerical 
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Figure  3.    Run  3,  Free  Surface  Deformation,  8  =  .09,  T*  ■  -1.1 
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Figure  5.    Run  5,  Free  Surface  Deformation,  5  =  .11,  T*  =  -1.05 
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instabilities.  Essentially,  any  increase  in  8  above  approximately  0.1 
serves  no  other  purpose  than  to  mask  a  number  of  surface  vortices 
unnecessarily,  with  only  a  trivial  gain  in  additional  time  steps  of 
execution.  Figures  8  through  14  represent  the  free-surface  deforma- 
tion at  execution-halt  time  and  clearly  demonstrate  that  the  gain  in 
program  execution  time  is  minimal.  The  numerical  model  is  surpris- 
ingly insensitive  to  this  parameter  once  it  has  accomplished  its  pur- 
pose of  removing  the  singularity.  The  qualitatively  chosen  optimum 
value  of  5  was  determined  to  be  0.1. 

TABLE  1 
VARIATION  OF  5  VERSUS  BREAKDOWN  TIME 


Breakdown 

Run  # 

5 

Time  (T*) 

1 

.05 

-1.4 

2 

.07 

-1.25 

3 

.09 

-1.1 

4 

.1 

-1.1 

5 

.11 

-1.05 

6 

.13 

-1.05 

7 

.15 

-1.0 

C      NUMBER  OF  FREE  SURFACE  VORTICES  (N) 
1.      Introduction 

The  number  of  free  surface  vortices  present  during  any  par- 
ticular run  is  automatically  calculated  by  a  subroutine  within  the 
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Figure  11.    Run  4,  Free  Surface  Deformation,  5 
Execution  Halt  T*  =  -1.00 


=  .10, 
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Figure  13.    Run  6,  Free  Surface  Deformation,  8 
Execution  Halt  T*  =  -0.95 


=  .13, 
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Figure  14.    Run  7,  Free  Surface  Deformation,  5 
Execution  Halt  T*  =  -0.95 


=  .15, 
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program  when  given  certain  desired  parameters.  The  operator 
directs  the  program  to  place  the  surface  vortices  after  deciding  on  the 
following  conditions: 

a   The  initial  distance  between  the  first  free  surface  vortex  and  the 
centerline. 

b.  The  position   (x/bo)  where  the  geometric  progressional  place- 
ment of  vortices  ends  and  constant  spacing  begins. 

c.  The  actual  distance  (x/bo)  between  vortices  after  the  point  in  (b) 
above  has  been  reached. 

The  details  of  the  calculation  of  vortex  spacing  is  not  particu- 
larly important  but  the  actual  decision  of  how  many  vortices  to  place 
along  the  free  surface  will  critically  affect  the  performance  of  the 
model.  As  discussed  earlier,  as  vortex  density  increases,  although 
leakage  decreases  (a  desirable  result),  the  velocity  impact  of  one  vor- 
tex upon  another  increases  and  rapidly  forces  the  breakdown  of  the 
model  prematurely. 
2.      Application 

Numerous  program  executions  were  completed  to  determine 
the  sensitivity  of  the  numerical  model  to  the  number  of  discrete  free 
surface  vortices.  Table  2  represents  the  results  of  a  portion  of  these 
executions. 

The  desingularization  parameter  (5)  was  adjusted  in  order  to 
maintain  the  criteria  established  earlier  (avoid  masking  vortices 
unnecessarily  yet  still  "cover"  the  breakdown  point).  This  adjustment 
of  5  was  necessary  to  ensure  that  the  model  remained  as  close  to  its 
natural  form  as  possible.     For  example,  if  5  was  maintained  at  0.1 


25 


during  run  #10,  it  would  effectively  mask  a  much  larger  number  of 
discrete  vortices  and  thus  introduce  an  unnecessary  artificiality  into 
the  execution  of  the  model.  Similarly,  if  5  was  maintained  at  0. 1 
during  run  #8,  the  singularity  in  the  velocity  expression  would  barely 
even  be  removed.  Figures  15  through  17  depict  the  formation  of  the 
free  surface  immediately  before  breakdown. 

TABLE  2 
VARIATION  OF  N  VERSUS  BREAKDOWN  TIME 


Breakdown 

Run 

5 

N 

Time  (T*) 

8 

0.18 

38 

-1.0 

9 

0.10 

76 

-1.1 

10 

0.085 

149 

-1.1 

The  results  of  these  executions  reveal  the  model's  insensitiv- 
ity  to  a  variable  number  of  discrete  free  surface  vortices.  The  only 
concern  remaining  to  the  operator  as  a  result  of  this  particular  study  is 
to  balance  the  leakage  versus  the  propensity  of  the  vortices  to  orbit 
about  each  other,  as  discussed  earlier,  and  to  consider  the  total  com- 
putational time  required  to  complete  each  run.  For  example,  Runs  #8 
and  #9  require  only  a  moderate  amount  of  time  to  complete,  but  Run 
#10  required  over  90  minutes  of  real-time  computation.  As  a  result, 
any  execution  with  N  >  149  must  be  executed  during  low-computer- 
use  time  only.    Qualitatively  weighing  these  factors,  it  was  determined 
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Figure  15.    Run  8  ,  Free  Surface  Deformation,  N  =  38 
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Figure  16.    Run  9,  Free  Surface  Deformation,  N  =  76 
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Figure  17.    Run  10,  Free  Surface  Deformation,  N  =  149 
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that  the  optimum  number  of  free  surface  vortices  was  76.  Average 
vortex  density  (N/xmax)  was  held  constant  throughout  the  remainder 
of  the  parametric  study. 

D.      MAXIMUM  X/B0  VALUE  (XMAX) 

1.  Introduction 

Very  closely  related  to  the  total  number  of  vortices  in  the 
model  is  the  total  length  of  free  surface  that  they  will  be  responsible 
for  representing.  It  is  not  desirable  to  have  the  vortices  represent  an 
extensive  stretch  of  free  surface  since  it  would  prolong  the  total  com- 
putation time  and  serve  no  meaningful  function.  It  is  also  undesirable 
to  restrict  the  extent  of  the  free  surface  so  severely  that  it  would  not 
truly  represent  nature  as  closely  as  possible.  The  simple  goal  of  the 
study  of  the  numerical  model's  sensitivity  to  a  changing  xmax  value  is 
to  confirm  the  theory  that  the  model  would  behave  similarly  at  differ- 
ent N,  provided  the  vortex  density  remained  relatively  constant. 

2.  Application 

Two  executions  were  completed  in  an  attempt  to  confirm  the 
theory  of  vortex  density  within  the  model.  Table  3  represents  the 
results  of  these  executions. 
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TABLE  3 
VARIATION  OF  XMAX  VERSUS  BREAKDOWN  TIME 


Run  # 

5 

N 

xmax 

Breakdown 
Time  (T*) 

11 
12 

.1 
.1 

76 

38 

10 
5 

-1.0 
-1.1 

From  Figures  18  and  19,  the  form  of  the  free  surface  demon- 
strated no  significant  change  when  xmax  was  decreased  from  10  x/b0 

to  5  x/b0-    Essentially,  there  is  no  major  advantage  in  using  a  value 

exceeding  about  5  x/b0  since  all  vortex  action  occurs  well  within  that 

value.    The  only  remaining  advantage  of  varying  this  parameter  is  to 

increase  or  decrease  real  CPU  time.    Obviously,  as  the  value  of  xmax 

increases,  computational  time  increases  as  well.    For  the  remainder  of 
this  parametric  study,  a  value  of  5  x/b0  was  utilized. 

E.      CENTER  OF  THE  SURFACE  VORTEX  DISTRIBUTION  (XV) 
1.      Introduction 

Elnitsky  [Ref.  10]  pointed  out  that  a  major  advantage  of  the 
Boundary  Integral  method  is  the  propensity  of  the  free  surface  vortices 
to  migrate  to  the  regions  where  the  radius  of  curvature  of  the  free 
surface  is  smallest.  The  surface  vortices  "roll  down"  the  hump  of  the 
deformed  free  surface  and  collect  in  the  region  of  the  scar.  Thus  one 
achieves  the  greatest  discretization  at  the  positions  along  the  free 
surface  where  it  is  most  needed. 
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Figure  18.   Run  11,  Free  Surface  Deformation,  zmaz  =  10 
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Figure  19.    Run  12,  Free  Surface  Deformation,  zmaz  =  5 
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As  discussed  earlier,  the  center  of  free  surface  vortex  distri- 
bution is  an  operator-selected  parameter.  The  algorithm  that  calcu- 
lates the  initial  position  of  each  free  surface  vortex  will  allow  the 
distribution  of  vortices  to  emanate  from  any  x/b0  position.  From  that 
selected  position,  a  subroutine  will  geometrically  place  vortices  to  the 
left  until  the  centerline  (x/b0  =  0)  has  been  reached  and  to  the  right 
until  the  x/b0  position  of  xmax  has  been  reached.  Since  the  vortex 
density  along  the  free  surface  is  greatest  at  the  center  of  distribution, 
by  moving  this  point  the  operator  can  affect  the  vortex  density  at  any 
point  along  the  free  surface  later  in  time.  Naturally,  as  one  progresses 
away  from  the  center  of  distribution,  this  effect  is  minimized  and  is 
eventually  eliminated  completely.  The  purpose  of  the  variation  of  this 
parameter  is  to  ascertain  whether  there  is  an  ideal  (xv)  position 
enabling  the  surface  vortices  to  "roll"  into  an  optimum  position  at  the 
correct  time  and  minimize  or  delay  the  onset  of  the  numerical  insta- 
bilities that  are  introduced  into  the  free  surface. 
2,     Application 

Runs  13  through  15  represent  a  sample  of  the  executions  of 
the  code  completed  in  order  to  identify  the  sensitivity  of  the  model  to 
the  parameter  (xv).  Table  4  lists  the  details  of  each  run  and  its  effect 
on  completion  time. 
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TABLE  4 
VARIATION  OF  xv  VERSUS  BREAKDOWN  TIME 

Breakdown 
Run#  xv(x/bo)  Time  (T*) 


13 

1.0 

-1.10 

14 

0.8 

-1.05 

15 

0.5 

-1.25 

Although  the  breakdown  time  of  runs  13  through  15  did  not 
vary  significantly  from  earlier  runs,  the  variation  of  xv  resulted  in 
totally  unacceptable  free  surface  formations.  From  Figures  20  through 
22,  it  is  noticed  that  the  flat  shape  of  the  hump  later  in  time  is  caused 
by  the  relatively  low  vortex  density  Also  noted  is  the  fact  that  the  free 
surface  vortices  continue  to  "roll  downhill"  but  terminate  farther  than 
in  previous  executions.  In  fact,  this  tendency  of  the  surface  vortices  to 
roll  farther  away  from  centerline  actually  contributes  to  the  formation 
of  the  numerical  instabilities  by  creating  sharp  corners  and  large 
spaces  between  individual  vortices.  Accordingly,  it  was  decided  to 
maintain  the  center  of  vortex  distribution  at  the  centerline. 
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Figure  20.    Run  13,  Free  Surface  Deformation,  xv  =  1.0 
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Figure  21.    Run  14,  Free  Surface  Deformation,  xv  =  0.8 
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Figure  22.   Run  15,  Free  Surface  Deformation,  xv  =  0.5 
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HI.   SMOOTHING 

A      INTRODUCTION 

Despite  all  efforts  to  minimize  the  effect,  or  delay  the  onset,  of 

the  numerical  Helmholtz  instabilities  during  the  parametric  study, 

they  consistently  occur  and  are  clearly  inevitable.     Extensive  efforts 

were  expended  to  discover  an  additional  parameter  within  the  model 

that  would   minimize  the   effect  of  the  instabilities  or  delay  their 

formation.     During  the  course  of  the  parametric  study,  as  a  path  of 

investigation,  it  was  decided  to  observe  the  variation  of  the  circulation 
strength  of  the  free  surface  vortices  with  respect  to  increasing  x/b0 

distance  as  measured  from  center  line.  It  was  theorized  that  sharp 
differences  in  circulation  values  between  adjacent  vortices  were 
responsible  for  the  tendency  of  the  adjacent  vortices  to  rotate  about 
each  another. 

As  a  representative  sample,  run  12  was  utilized  to  depict  the  vari- 
ation of  the  circulation  with  respect  to  x/b0.  Figure  23  clearly  demon- 
strates the  so  called  "sawtooth"  instabilities.  Each  point  on  the  figure 
represents  an  individual  free  surface  vortex  circulation  strength. 

It  is  useful  to  observe  the  interaction  of  the  Helmholtz  instabilities 
with  the  variation  of  the  circulation  values.  Figure  24  depicts  the  free 
surface  vortices  that  eventually  (within  one  or  two  time  steps)  will  lead 
to  the  onset  of  a  Helmholtz  instability.     Figure  25  represents  the 
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Figure  23.    Run  12,  Variation  of  V  With  Respect  to  x/bQ 
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Figure  24.  Run  12,  Discrete  Vortices  Leading  to  Instabilities 
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Figure  25.    Run  12,  Variation  of  V  With  Respect  to  x/bQ 
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corresponding  value  of  circulation  for  each  free  surface  vortex  at  the 
time  associated  with  Figure  24.  The  sharp  variation  in  circulation 
between  the  vortices  is  noted.  At  the  actual  point  of  breakdown,  the 
difference  in  circulation  effectively  approaches  infinite  values.  The 
purpose  of  attempting  to  smooth  the  variation  of  circulation  is  to  delay 
the  inevitable  formation  of  the  Helmholtz  instability. 

B.     APPLICATION 

From  Longuet-Higgins  [Ref  12],  a  five-point  smoothing  routine  was 
developed  and  implemented  to  smooth  the  value  of  the  circulation 
when  a  set  of  adjustable  parameters  have  been  established  and  a  set  of 
specific  conditions  have  been  met.  With  a  minor  modification  of  the 
code,  the  operator  can  elect  to  either  activate  the  smoothing  subrou- 
tine continuously  or  require  that  a  specified  condition  be  met  prior  to 
commencing  the  smoothing  routine.  Additionally,  the  operator  may 
specify  the  exact  time  where  the  routine  commences  looking  for  the 
desired  conditions.  In  this  particular  investigation,  the  smoothing 
routine  was  called  when  the  difference  in  circulation  values  between 
two  adjacent  vortices  was  greater  than  0.3  and  the  nondimensional 
time  was  at  least  -2.1. 

Other  variables  that  the  operator  can  adjust  include: 

1.  Start  time.  The  operator  has  the  ability  to  specifically  state 
when,  in  nondimensional  time,  the  smoothing  routine  is  called. 
Clearly  it  serves  no  useful  purpose  to  smooth  at  early  times 
because  the  free  surface  has  not  started  deforming.  Since 
smoothing  is  an  undesirable  solution  of  the  problem  in  any  event, 
its  use  should  be  avoided  for  as  long  as  possible.  A  nondimen- 
sional time  of  -2.1  was  chosen  as  providing  a  comfortable  margin 
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of  smoothing  yet  not  unnecessarily  introducing  artificiality  into 
the  problem. 

2.  Frequency  of  smoothing.  Once  the  actual  start  time  of  smooth- 
ing has  been  reached,  the  operator  may  direct  how  often  future 
smoothing  attempts  will  be  made.  Arbitrarily,  it  was  decided  that 
once  smoothing  started,  it  would  occur  every  25  time  steps. 

3.  Duration  of  smoothing.  Once  smoothing  has  occurred  for  a  spe- 
cific time  step,  the  operator  may  (with  a  minor  modification  to 
the  code)  elect  to  automatically  smooth  the  next  "N"  time  steps 
prior  to  starting  the  smoothing  counter  once  again.  This  selec- 
tion is  useful  in  that  the  Longuet-Higgins  method  is  not  a  domi- 
nating technique,  and  the  circulation  variation  may  be  so  large 
that  it  requires  several  consecutive  smoothing  time  steps  to 
remove  the  sharp  difference  which  initiated  the  routine.  In  this 
application  only,  "N"  represents  the  number  of  consecutive  time 
steps  that  smoothing  will  be  completed. 

4.  X/b0  position.  It  is  unnecessary  to  smooth  the  variation  of  the 
circulation  at  points  along  the  free  surface  where  instabilities  do 
not  exist  or  normally  form.  In  general,  the  instabilities  form  at  an 
x/b0  range  of  (0.5  <  x/b0  <  1.25)  at  or  very  close  to  the  formation 
of  the  scar. 

Figure  26  represents  a  duplication  of  run  12  depicted  in  Figure 

23  except  with  the  Longuet-Higgins  continuous  smoothing  routine 

activated.     For  this  particular  run,  the  following  smoothing  routine 

selections  were  made: 

start  time:  -2.1 

smoothing  duration:  five  time  steps 

time  delay  25  time  steps 

x/bo  range  0.5  <  x/b0  <  1.25 

Notice  that  the  sharp  peak  has  been  reduced  yet  the  basic  form  of  the 
free  surface  deformation  has  been  retained. 
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Figure  26.    Run  12,  Smoothed  Variation  of  V  With  Respect  to  x/bQ 
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IV.    RE-DISCRETIZATION 

A      INTRODUCTION 

After  close  scrutiny  of  each  run  during  this  parametric  study,  it 
was  determined  that  developing  a  method  of  re-discretizing  the  free 
surface  would  be  useful.  This  task  required  the  development  of  a 
routine  within  the  model  that  would  add  a  discrete  vortex  at  a  specific 
location.  The  critical  problem  would  be  to  determine  exactly  when 
and  where  to  add  a  free  surface  vortex  element  and  what 
corresponding  circulation  to  assign  to  it.  The  routine  must  have  the 
ability  to  interrupt  program  execution  during  any  specific  time  step, 
make  the  determination  that  a  new  vortex  is  required,  and  place  it 
while  still  maintaining  the  continuity  of  the  program  and  the  general 
form  of  the  free  surface. 

R      APPLICATION 

A  subroutine  was  developed  within  the  code  and  assigned  an 
activation  option.  If  the  "add  vortex"  option  is  selected,  the  program 
will  add  a  discrete  free  surface  vortex  element  under  specific 
conditions.  This  routine  calculates  the  complex  distances  between 
each  free  surface  vortex  element  and  stores  the  information  in  an 
array.  When  the  difference  of  the  rate  of  change  of  distance  between 
two  adjacent  vortices  exceeded  an  adjustable  percentage  value,  the 
subroutine  "ADDVOR"  was  called.  For  the  purposes  of  this 
investigation,  the  specific  value  utilized  was  33  percent. 
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An  additional  restriction  was  imposed  upon  the  routine  requiring 
that  the  two  vortices  under  consideration  be  increasing  the  distance 
between  them  instead  of  closing  the  distance  between  them.  If 
indeed  the  distance  between  two  surface  vortices  was  closing  at  an 
acceptable  rate  (i.e.,  greater  than  the  33  percent  requirement),  the 
subroutine  would  add  vortices  at  those  positions  along  the  free  surface. 
This  is  obviously  undesirable  because  the  vortices  are  already  closely 
packed.  Figure  27  represents  a  typical  situation  where  a  vortex  would 
be  added. 

Operator  considerations  include: 

1.  Closing  rate.  As  discussed  earlier,  the  rate  of  change  of  distance 
between  vortices  will  dictate  when  the  program  will  add  a  vortex. 
There  is  no  restriction  as  to  the  percentage  that  is  chosen  but  it 
was  determined  that  an  acceptable  range  was  between  30 
percent  and  50  percent. 

2.  Circulation.  The  value  of  the  circulation  of  the  added  vortex  was 
calculated  algebraically  between  the  two  vortices  so  that  the  slope 
of  the  circulation  with  respect  to  x/bQ  would  be  preserved 
between  the  (now)  three  surface  vortices.  The  values  of  the  new 
circulations  are  given  by  the  following  expressions: 

r;  .(ST.-r.j/e  (5) 


^2     =<5*r2    "ri)/6 


(6) 


rncw  =  (r1  +  r2)/3  (7) 
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Figure  27.   Typical  Free  Surface  Condition  Requiring 
an  Additional  Vortex 
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Tl  represents  the  circulation  of  vortex  #1  and  r2  represents  the 
circulation  of  vortex  #2  prior  to  the  addition  of  the  new  vortex. 
r,  and  T2  represent  the  new  calculated  vortex  circulation 
values,  respectively.  rnew  represents  the  value  of  the  circulation 
of  the  new  vortex. 

3.  New  vortex  position.  In  order  to  accurately  place  the  new  vortex, 
it  was  decided  to  calculate  the  new  Z  position  utilizing  a 
conservation  of  circulation  principle.  Accordingly,  the  position  of 
the  new  vortex  (Znew)  is  given  by  the  expression: 


znew  =  ((r1  *z1  +  r2  *z2)-(r;  *Zl  +  r;  *z2))/r, 


(8) 


Although  several  placement  and  circulation  algorithms  were 
developed  and  tested,  the  ones  listed  in  equations  (5)  through 
(8)  performed  the  most  adequately.  Figure  28  represents  the 
execution  of  this  routine  immediately  following  the  addition  of  a 
vortex.  This  figure  follows  Figure  27  in  time.  Notice  that  the 
routine  placed  the  vortex  in  a  satisfactory  manner.  Figure  29 
represents  the  circulation  strength  assigned  to  the  added  vortex. 

Smoothing  As  can  be  seen  from  Figure  29,  although  the 
circulation  is  correctly  assigned  to  the  new  vortex,  there  is  a 
slight  difference  between  the  new  value  and  the  circulation  values 
of  its  adjacent  vortices.  Accordingly,  the  operator  may  elect  to 
immediately  smooth  the  circulation  values  after  adding  a  vortex. 
This  may  not  be  necessary,  however,  if  the  value  of  the  circulation 
difference  exceeds  0.3  of  its  neighbors  as  mentioned  earlier. 
Figure  30  depicts  a  smoothed  version  of  the  circulation  strength 
of  the  added  vortex. 
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Figure  28.    Free  Surface  Deformation  Immediately  After 

Adding  a  Vortex 
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Figure  29.    Circulation  Value  of  an  Added  Vortex 
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Figure  30.    Smoothed  Circulation  Value  of  an  Added  Vortex 
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V.  STRATIFIED  FLOW  EXPERIMENTS 

A   EXPERIMENTAL  APPARATUS 

Specific  details  of  the  experimental  apparatus  were  discussed 
extensively  by  Leeker  [Ref  11]. 

Experiments  were  conducted  in  a  water  basin  approximately 
twelve  feet  long,  three  feet  wide,  and  four  feet  deep.  Two  counter- 
rotating  vortices  were  generated  through  the  use  of  two  opposing 
plates  activated  by  a  system  of  weights  and  pulleys.  Flow  visualization 
was  performed  through  the  use  of  a  fluorescene  dye  injected  through 
openings  on  the  underside  of  each  plate.  A  video  camera,  a  VHS  video 
cassette  recorder,  a  monitor,  and  a  date/timer  were  utilized  as  sup- 
plementary equipment. 

R   PROCEDURES 

Basic  experimental  procedures  were  thoroughly  discussed  by 
Leeker  [Ref  11].  The  procedures  in  this  investigation,  although  fun- 
damentally similar,  differed  from  those  of  Leeker  with  respect  to  the 
addition  of  a  sharp  salt-water/fresh-water  density  interface. 

To  obtain  a  sufficiently  sharp  density  difference,  50  pounds  of 

ordinary  salt  were  thoroughly  mixed  with  465.7  gallons  of  fresh  water 

in  the  tank  bottom  and  allowed  to  settle  overnight. 

From  Striftos  [Ref  14],  the  concentration  is  given  by: 

weight  [water] 


1/C=1  + 


weight  [salt]  jgj 
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The  density  of  the  salt  water  layer  is  given  by  the  expression: 

Ps=Po*U+.71*C)  (10) 

Where  p0  is  the  density  of  the  fresh  water  (62.43  lbs/ft3). 

The  density  difference  then  is  simply  given  by  (Ps-po).  For  the 
purposes  of  these  experiments,  C  =  .0127  lbs.  salt/lbs.  water,  and 
ps  =  62.993  lbs./ft3.   Finally,  Ap  =.563  lbs./ft3. 

A  freshwater  dye  combination  was  introduced  by  gravity  feed  from 
a  reservoir  through  small  tubes  which  spilled  onto  a  floating  board  in 
the  tank.  The  dyed  freshwater  trickled  over  the  edge  of  the  board  and 
settled  on  top  of  the  salt-water  layer.  The  tank  was  filled  in  this  man- 
ner for  several  hours,  resulting  in  a  sharp  density  interface.  The 
remainder  of  the  experiments  were  carried  out  in  the  normal  method. 
Efforts  were  concentrated  on  recording  the  interaction  of  the  vortex 
pair  with  the  interface. 

C   RESULTS  OF  THE  EXPERIMENTAL  RUNS 

Approximately  ten  experimental  runs  were  completed.  Two  runs 
(04  02  and  04  08)  were  chosen  to  discuss  as  representative  samples. 
Analysis  procedures  of  all  experimental  runs  were  identical  to  those 
discussed  by  Leeker  [Ref  11]. 

Figure  31  depicts  the  right  vortex  (and  image)  center's  path 
for  run  04  02.  Figure  32  depicts  the  nondimensional  velocity  profile 
of  the  right  vortex  center  also  for  run  04  02.  Based  on  the  maximum 
velocities  displayed  by  Figure  32,  the  Froude  number  for  the  vortex  in 
saltwater  =  0.78  and  for  freshwater  =  1.05. 
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Figure  31.    Vortex  Center  Path  for  Run  04  02,  Fr  =  1.05 


55 


1.2 


1.0    . 


0.8    ■ 


0.6 


0.4    . 


0.2 


0.0 


Figure  32.   Vortex  Center  Nondimensional  Velocity 
for  Run  04  02,  Fr  =  1.05 
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From  Figure  31,  the  vortex  center  clearly  "saw"  the  density 
interface  as  a  makeshift  free  surface  and  began  to  deflect.  Upon  intru- 
sion into  the  lower-density  fluid,  the  velocity  increased  and  the  vortex 
path  became  more  vertical;  in  effect,  the  vortex  no  longer  "saw"  a 
barrier.  As  the  vortex  center  is  influenced  by  the  actual  free  surface,  it 
behaved  much  like  previous  investigations  [Refs  10,  11].  The  two  dis- 
tinctly different  maximum  velocity  peaks  on  Figure  32  support  this 
observation. 

The  vortex  in  run  04  08  was  only  visible  while  in  the  lower 
salt-water  layer.  Figure  33  depicts  the  vortex  and  image  path.  Figure 
34  represents  the  nondimensional  velocity  profile  of  the  vortex  cen- 
ter. In  this  run,  although  the  center  was  not  tracked  through  the 
interface,  the  vortex  also  appeared  to  "see"  that  there  was  a  slight  free 
surface  above  it,  as  evidenced  by  the  outward  slant  of  its  path. 
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Figure  33.    Vortex  Center  Path  for  Run  04  08,  Fr  =  .38 
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Figure  34.   Vortex  Center  Nondimensional  Velocity 
for  Run  04  08,  Fr  =  .38 
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VI.   CONCLUSIONS 

The  numerical  analysis  has  shown  that  the  evolution  of  the  free 
surface  can  be  calculated  through  the  use  of  the  line  vortices.  No 
filtering  is  necessary  provided  that  the  complex  velocity  function  is 
desingularized  in  a  suitable  manner.  The  analysis  becomes  more 
robust  as  the  Froude  number  is  increased.  Even  then,  the  inception  of 
small-wavelength  Helmholtz  instability  cannot  be  avoided.  However, 
the  time  at  which  such  instabilities  occur  can  be  delayed  by  a  suitable 
arrangement  of  the  initial  position  of  the  vortices.  Fortunately, 
however,  the  critical  time  at  which  the  instability  manifests  itself  is 
near  the  time  at  which  the  wave  breaking  is  likely  to  occur.  This 
analysis  is  not  equipped  to  deal  with  wave  breaking.  For  smaller 
Froude  numbers  (say  Fr  -  0.50),  the  agreement  between  the  calculated 
and  measured  free  surface  shapes  is  quite  encouraging.  A  better 
agreement  would  have  been  a  source  of  puzzlement  rather  than  a 
source  of  gratification  because  the  vortices  in  the  experiments  become 
quickly  turbulent,  diffuse  over  a  large  area  and  lose  their  strength  due 
to  the  cancellation  of  oppositely  signed  vorticity  in  the  overlapping 
regions  of  vortices  and  due  to  the  negative  vorticity  generated  at  the 
free  surface.  Although  the  numerical  model  has  the  capability  to  vary 
the  vortex  strength  with  time  (e.g.,  using  a  diffusing  Oseen  vortex), 
this  was  not  done. 
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For  any  of  the  Froude  numbers  encountered  in  the  analysis  or 
experiments,  no  wave  train  was  observed  on  either  side  of  the  scars. 
This  conclusion  was  somewhat  anticipated  on  the  basis  of  previous 
calculations.  For  very  small  Froude  numbers  (say,  Fr  <  0.15),  the  scars 
were  very  small,  the  vortices  followed  the  rigid-wall  path,  and  the 
scars  were  slaved  to  the  vortices,  as  found  previously  by  Sarpkaya  and 
Henderson  [9]. 

The  numerical  model  can  be  improved,  within  the  limitations  of 
the  inviscid  flow  theory,  by  using  vortex  sheet  panels  rather  than  line 
vortices.  It  is  also  possible  to  use  an  adaptive  grid  system  together 
with  a  suitable  finite- difference  method  to  solve  the  equations  of 
motion  directly.  Nevertheless,  the  difficulties  stemming  from  the 
occurrence  of  Helmholtz  instability  during  the  roll-up  of  the  vortex 
sheets,  axial  instability  of  the  vortex  pair,  three-dimensional 
instabilities  on  the  free  surface  and  the  turbulent  diffusion  of  the 
vortices  remain  as  the  most  serious  impediments  towards  an  exact 
solution  of  the  problem. 

In  summary,  the  principal  objectives  of  this  investigation  have 
been  accomplished.  The  extensive  parametric  study  has  shown  that 
the  numerical  model  is  relatively  insensitive  to  a  variety  of  changing 
parameters  and  that  the  Helmholtz  instabilities  are  inevitable  in  the 
execution  of  this  code  beyond  some  critical  time. 

The  supplementary  objective  of  analyzing  the  interaction  of  a 
vortex  pair  with  a  sharp  density  interface  was  also  accomplished. 
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Although  only  a  preliminary  to  major  investigations  in  the  future,  this 
"first  look"  of  the  interaction  of  two  counter-rotating  vortices  with  a 
sharp  density  interface  revealed  interesting  information,  such  as  the 
considerable  increase  in  velocity  of  the  vortex  pair  after  its  intrusion 
into  the  fresh-water  layer. 
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